Fitting and Choosing Models
\[y(t) = m(t) + e(t)\]
\(m(t)\) is the deterministic component, which we can explain by, e.g., a constant increase (linear trend), a periodic function, a combination of both, …
\(e(t)\) is the stochastic component, that may or may not be random (white noise)
Assuming this is a sinus function, \[\alpha_1 + \alpha_2 \sin(t + \alpha_3)\], we need non-linear regression (function is non-linear in \(\alpha_3\))
A possible approach is to find the minimum for Akaike’s Information Criterion (AIC) for ARMA(\(p,q\)) models and series of length \(n\): \[AIC = \log \hat{\sigma}^2 + 2(p+q+1)/n\] with \(\hat{\sigma}^2\) the estimated residual (noise) variance.
Instead of finding a single best model using this single criterion, it may be better to select a small group of “best” models, and look at model diagnostics for each: is the residual white noise? does it have stationary variance?
Even better may be to keep a number of “fit” models and consider each as (equally?) suitable candidates.
[1] -23547.93 -30235.42 -30713.51 -30772.31 -30815.14 -30816.35 -30818.27
[8] -30818.39 -30817.82 -30815.84
[1] -23994.90 -30320.26 -30843.40 -30885.36 -30945.22 -30943.51 -30952.06
[8] -30957.84 -30955.86 -30955.38
# Prediction, e.g. with AR(6)
x = arima(T.outside, c(6,0,0))
pltpred = function(xlim, Title) {
plot(T.outside, xlim = xlim, type='l')
start = nrow(meteo) + 1
pr = predict(x, n.ahead = xlim[2] - start + 1)
lines(start:xlim[2], pr$pred, col='red')
lines(start:xlim[2], pr$pred+2*pr$se, col='green')
lines(start:xlim[2], pr$pred-2*pr$se, col='green')
title(Title)
}
pltpred(c(9860, 9900), "10 minutes")# Prediction, e.g. with AR(6)
T.per = 18.2-4.9*sin(pi*(hours+1.6)/12)
an = T.outside - T.per
x2 = arima(an, c(6,0,0))
pltpred2 = function(xlim, Title) {
plot(an, xlim = xlim, type='l')
start = nrow(meteo) + 1
pr = predict(x2, n.ahead = xlim[2] - start + 1)
lines(start:xlim[2], pr$pred, col='red')
lines(start:xlim[2], pr$pred+2*pr$se, col='green')
lines(start:xlim[2], pr$pred-2*pr$se, col='green')
title(Title)
}
pltpred2(c(9860, 9900), "10 minutes")Now that we have an idea of the trend and the anomaly, we can combine them to make predictions
plot(T.outside,xlim=c(1,19970), type='l')
x.an = arima(an, c(6,0,0)) # model the anomaly by AR(6)
x.pr = as.numeric(predict(x.an, 10080)$pred)
x.se = as.numeric(predict(x.an, 10080)$se)
hours.all = c(meteo$hours, max(meteo$hours) + (1:10080)/60)
T.per = 18.2-4.9*sin(pi*(hours.all+1.6)/12)
lines(T.per, col = 'blue')
hours.pr = c(max(meteo$hours) + (1:10080)/60)
T.pr = 18.2-4.9*sin(pi*(hours.pr+1.6)/12)
lines(9891:19970, T.pr+x.pr, col='red')
lines(9891:19970, T.pr+x.pr+2*x.se, col='green')
lines(9891:19970, T.pr+x.pr-2*x.se, col='green')
title("predicting 1 week: periodic trend + ar(6) for anomaly")We can also use these models to simulate data
Prediction/forecasting:
A phase shift model \(\alpha \sin(x + \phi)\) can also be modelled by
\(\alpha sin(x) + \beta cos(x)\); this is essentially a re-parameterization.
A phase shift model \(\alpha \sin(x + \phi)\) can also be modelled by
\(\alpha sin(x) + \beta cos(x)\); this is essentially a re-parameterization.
Try the following code:
x = seq(0, 4*pi, length.out = 200) # x plotting range
sincos = function(dir, x) { # plot the combination of a sin+cos function, based on dir
a = sin(dir)
b = cos(dir)
# three curves:
plot(a * sin(x) + b * cos(x) ~ x, type = 'l', asp=1, col = 'green')
lines(x, a * sin(x), col = 'red')
lines(x, b * cos(x), col = 'blue')
# legend:
lines(c(10, 10+a), c(2,2), col = 'red')
lines(c(10, 10), c(2,2+b), col = 'blue')
arrows(x0 = 10, x1 = 10+a, y0 = 2, y1 = 2+b, .1, col = 'green')
title("red: sin(x), blue: cos(x), green: sum")
}
for (d in seq(0, 2*pi, length.out = 100)) {
sincos(d, x)
Sys.sleep(.1)
}So, we can fit the same model by a different parameterization:
…which is a linear model, identical to:
nlm()The EVI (enhanced vegetation index) is a one-dimensional index computed from near infrared, red, and blue reflection measurements. Values can be interpreted as “greenness” and range from -1 to 1.
For each of the \(1200^2\) pixels, we interpret the set of observations as time series.
Consider the general model: \[y_t = T_t + S_t + R_t\]
We could be interested in both, changes in the trend and seasonality.
An example series (1 pixel in forest area):
We compute differences from the overall mean to center the series around 0.
We compute differences from the overall mean to center the series around 0.
We assume yearly seasonality and build seasonal sub-series:
\[ \begin{aligned} y_t^{(16)} &= (y_{2005,16},y_{2006,16},y_{2007,16},y_{2008,16}, ...) \\ y_t^{(32)} &= (y_{2005,32},y_{2006,32},y_{2007,32},y_{2008,32}, ...) \\ ... \end{aligned} \] for which we simply compute mean values.
After subtracting the seasonal component, the series looks like:
After removing large parts of seasonality, we can study the trend of the data. Deforestation should be visible by abrupt changes in the trend structure. To test, whether there is a structural change, we assume the following linear regression model \[T_t = \beta_0 + \beta_1 t\] and hypothesize that the regression coefficients \(\beta_0, \beta_1\) do not change in time.
Fitting a global trend leads to:
Looking at the cumulated sum of the model’s residuals shows that they do not fluctuate around 0, which we expected in our hypothesis.
By thresholding the cumulated residuals based on a given confidence level, we can find a breakpoint around t=107.
A better model:
Verbesselt, J., Hyndman, R., Newnham, G., & Culvenor, D. (2010). Detecting trend and seasonal changes in satellite image time series. Remote Sensing of Environment, 114, 106-115. DOI: 10.1016/j.rse.2009.08.014.
R. B. Cleveland, W. S. Cleveland, J.E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Journal of Official Statistics, 6, 3-73.
W. S. Cleveland, E. Grosse and W. M. Shyu (1992) Local regression models. Chapter 8 of Statistical Models in S eds J.M. Chambers and T.J. Hastie, Wadsworth and Brooks/Cole.